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I. INTRODUCTION 



There has been intense research interest in phase transitions in mass-transport and growth mod- 
els involving adsorption and desorption, fragmentation, diffusion and aggregation. These pro- 
cesses are ubiquitous in nature and arise in a large number of seemingly diverse systems such as 
growing interfaces Jl[Q], colloidal suspensions |3], polymer gels [4], river networks [5], granular 
materials Q6Q, traffic flows [|70, etc. In these systems, different nonequilibrium states arise if the 
rates of the underlying microscopic processes are varied. Conservation laws also play an important 
role in determining both the time-dependent behavior and steady states of such systems. As these 
steady states are usually not described by the Gibbs distribution, they are hard to determine. How- 
ever, much insight on this issue has been gained by studying lattice models. Due to their simplistic 
nature, these models can be treated either exactly or via a mean-field (MF) approach [|8|4l3n. They 
are also simple to implement numerically using Monte Carlo (MC) techniques [|14l ll5ll. 



In this paper, we present a pedagogical discussion of the modeling and simulation of mass- 
transport and growth phenomena. We discuss analytical and numerical techniques in the context 
of mass-transport models where the elementary move is the fragmentation of mass k, and its 
subsequent diffusion to a neighboring site where it aggregates. The A;-chip models that we study 
here are interesting in physical situations where the deposited material consists of polymers. We 
study the MF limit of these models, focusing on the steady-state mass distribution [P(m) vs. m], 
which is characterized by k branches. We also compare the MF results with MC simulations in d 
= 1,2. 

This paper is organized as follows. In Sec. UH we present a framework for mass-transport 
models in terms of the rate (evolution) equations for P(m, t), the probability that a site has mass 
m at time t. We then discuss various systems which can be described within this framework. In 
Sec. Unl we introduce A;-chip models and obtain analytical results for the MF versions of these 
models. In Sec. [TV] we present MC results for mass-transport models, and compare them with the 
corresponding MF solutions. We conclude this paper with a summary and discussion in Sec. |V] 
The appendices contain details of calculations and MC procedures. 
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II. FRAMEWORK AND APPLICATIONS OF MASS-TRANSPORT MODELS 



A. Lattice Models and Rate Equations 

We consider lattice models of mass transport with the processes of fragmentation, diffusion 
and aggregation. For simplicity, we describe the models on a 1 -dimensional lattice with periodic 
boundary conditions. (The generalization to higher dimensions is straightforward.) To begin with, 
masses are placed randomly at each site with an overall mass density p. Let mj(i) be the mass at 
site i at time t. The mass variables assume discrete values 0,1,2,3, etc. The evolution of the system 
is as follows. A piece of mass n chips off a site having mass m (> n) with rate g m (n). This piece 
deposits on the right neighbor with probability p, or on the left neighbor with probability 1 — p. 
The mass of the chosen neighbor adds up, while that of the departure site decreases, with the total 
mass of the system remaining conserved. Figure 1 is a schematic depiction of the above model. 
To facilitate MC simulations, the update rules can be rewritten as follows: 

1) Randomly pick a site i at time t with mass rriiit) = m. The site is updated as rrii(t + 1) = 
rriiit) — n with rate g m (n). 

2) The neighboring sites are updated as m i+ i [t + 1) = m i+1 (t) + n with probability p, or m.j_i (t + 
1) = mi-i{t) + n with probability 1 — p. 

We also study the above models within a MF approximation which keeps track of the distribu- 
tion of masses, ignoring correlations in the occupancy of adjacent sites. Although the MF theory 
has this shortcoming, our MC simulations show that it gives an accurate description of the above 
model, even in the 1 -dimensional case. Let P(m, t) denote the probability that a site has mass m 
at time t. In the MF limit, P(m, t) evolves as follows: 

m oo m,2 

Pirn. t\ = -P(m t\ n.„Jm.A - P(m t\ P(m.„ t\ 

dt 



—P(m,t) = -P(m,t) 9m(mi) - P(m,t) ^ p ( m 2,t) ^ 9 m2 ( m i) 

mi=l rri2=l mi=l 



+ 2j P(m + m 1 ,t)g m+mi (m 1 ) 



mi=l 



+ y ] P{m - m^t) P(m 2 ,t)g m2 (m 1 ), m > 1, (1) 

mi=l m,2=n\\ 



d °° 

-P(0,t) = -P{0,t) P(m 2 ,t)J29m 2 (mt)+ Yl P(mi,t)g mi ( mi ). (2) 

7712=1 mi=l mi=l 



7712 

p(n t\ = -pen t\ pc m „ t\ 

dt 

These equations enumerate all possible ways in which a site with mass m may change its mass. 
The first term on the right-hand- side (RHS) of Eq. (OQ) is the "loss" of mass m due to chipping, i.e., 



a site with mass m may lose a fragment of mass mi (< m) to a neighbor. The second term on the 
RHS represents the loss due to transfer of mass from a neighbor chipping. The third and fourth 
terms are the "gain" terms which represent the ways in which a site with mass greater (lesser) than 
m can lose (gain) the excess (deficit) to yield mass m. The terms of Eq. © can be interpreted 
similarly. In order to ensure that all loss and gain terms have been included in the rate equations, 
it is useful to check the sum rule 

d °° 00 

-^P(m,t) = 0, or ^P(m,t) = l. (3) 

m=0 m=0 

With some algebra, it can be shown that Eqs. ©-© do indeed satisfy the above rule. We provide 
the steps for this check in Appendix A. 

If other microscopic processes are present, additional terms will have to be included in the rate 
equations (OQ) and ©. For example, if adsorption of a unit mass at a site occurs with rate q, we 
require additional terms —qP(m, t) and +qP(m— 1, t) in Eq. ©, and —qP(0, t) in Eq. ©. While 
it is simple to write realistic rate equations by including all relevant microscopic processes, it is 
often difficult to solve them analytically to obtain either time-dependent or steady-state solutions. 

Several MF models studied earlier may be obtained as special cases of Eqs. ©-© by an ap- 
propriate choice of the chipping kernel g m (n). Some interesting issues which have been addressed 
in these studies, in addition to obtaining steady-state mass distributions, are the possibilities of 
phase transitions in these models. We mention two representative examples to highlight the typi- 
cal questions which are addressed in this area. 



1) Majumdar et al. Ill ill studied a conserved-mass model in which either a single unit or the entire 



mass could dissociate from a site. Thus, g m ( n ) has the form 

g m (n) = w5 nA + 5 n>m , (4) 

where w is the relative rate of the 1-chip process. The corresponding rate equations are obtained 
by substituting Eq. © in Eqs. CD)-© as follows: 

—P(m,t) = -(l+w)(l + s 1 )P(m,t)+wP{m + l,t)+wsiP(m-l,t) 

m 

+ P(m- m 1 ,t)P(m 1 ,t), rn > 1, (5) 

mi=l 

j t P(0, t) = -(1 + w) Sl P(0, t) + wP(l, t) + Sa , (6) 
where s 1 = YZ=\ P ( m > *)• 



The steady-state mass distributions [P(m) vs. m] for Eqs. ©-(|6]) were calculated by Majum- 
dar et al. as a function of the density p = (to) of the system. The relevant analytical techniques 
are described in Sec. [Till They observed a dynamical phase transition as p was varied (w being 
fixed), with the different phases being characterized by different steady-state distributions. For 
p < p c (w), P(m) decayed exponentially for large m. For p = p c (w), P(m) showed a power-law 
decay, P(m) ~ vrT T with a universal exponent r = 5/2. Finally, the "high-density" phase arising 
for p > p c (w) was characterized by the formation of an infinite aggregate (at m = oo). The 
aggregate coexisted with smaller clusters, and their mass distribution showed a power-law decay, 
P(m) ~ mT T . 

2) Rajesh et al. J12I1 studied a system of fragmenting and coagulating particles with mass- 
dependent diffusion rates. In this model, 

g m (n) = w5 n ,i + m~ a 5 ntm . (7) 

The case a = corresponds to the model in Eq. ©. The corresponding rate equations are 

^-P(m, t) = — (w + m~ a + ws\ + Si) P(m, t) + wP(m + 1, t) + ws 1 P(m — 1, t) 
dt 



P(m — TOi, t)P(mi, t) 
ml 

: P{0,t) = -(ws 1 + s 1 )P(0,t) + wP(l,t) + s 1 , (9) 



+ ynm-m 1 ,t )nmi ,t) 

^— ' TO? 

mi=l 1 

d 



dt 

where s x = Ylm=i m ~ a ^ > ( m ^)- 

For a > 0, Rajesh et al. showed that there is no dynamical phase transition. The high-density 
phase with an infinite aggregate disappears, although its imprint in the form of a large aggre- 
gate is observed in finite systems. Further, the steady-state probability distribution P(m) decays 
exponentially with to for all p and a > 0. 

Before concluding this discussion, a few words regarding the condensation transition are in 
order. The condensation observed in the model in Eqs. ©-© occurs due to the dynamical rules 
of evolution and not due to an "attraction" between the masses. Though it shares analogies with 
Bose-Einstein condensation (BEC), an important difference is that these condensates are formed 
in real space rather than momentum space as is the case in BEC. As a matter of fact, condensation 
occurs in a variety of seemingly diverse systems which are governed by nonequilibrium dynamics 
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B. Some Applications of Mass-Transport Models 



We now discuss some physical applications of mass-transport models. The aim here is to stress 
the general nature of the questions addressed in a variety of physical situations involving mass 
transport. 



1. Magnetic Nanoparticles 

Recently, there has been much research interest in suspensions of single-domain magnetic 
nanoparticles (MNP), which have a wide range of technological applications, e.g., memory de- 



17] 



vices, magnetic resonance imaging, targeted drug delivery, bio-markers and bio-sensors R16 . 
A major reason for the utility of MNPs is the ease with which they can be detected and manipulated 
by an external magnetic field. Their response times are strongly size-dependent, thus introducing 
the possibility of controlling particle sizes to obtain desired response times. 

An inherent property of MNP suspensions is cluster formation, due to the presence of attractive 



interactions of varying strengths between the constituent particles [[18|]. Therefore, mass-transport 
models with fragmentation and aggregation have been traditionally employed to study cluster- 
ing dynamics in these systems. The steady-state cluster-size distributions and the average cluster 
size are determined by the interplay between aggregation (due to attractive interactions) and frag- 
mentation (due to repulsive interactions and thermal noise) [18]. Assuming that the number of 
particles is N, and denoting the number of clusters containing k particles at time t by c(k, t), the 



rate equations in the MF approximation are as follows [|19jl 



d 1 -v ^ 

_c(A;, t) = - S k ,i+jKijc(i, t)c(j, t) - c(k, t) ^ K kjc(j, t) 

i,j=l j=l 



oo 



+f k+ ic{k + 1, t) - f k c(k, t) + 8 k>1 fjc{j, *)> k>l. (10) 

In Eq. (TTOl) . Kij and f k are the aggregation and fragmentation kernels, respectively. The aggrega- 
tion kernel describes the coalescence of two clusters containing i and j particles to yield a larger 
cluster with k = i +j particles. In many models, it is assumed to have a mass-dependent form, 
= D^^+j^^). This accounts for the reduced mobility of large clusters. The fragmentation kernel 
fk describes the loss of one particle from a cluster with k particles, and has also been assumed to 
have a mass-dependent form, f k = wk u . Equation (TTOl) can be rewritten in terms of probability dis- 
tributions [cf. Eqs. CO)-©] by introducing a normalization factor, P(k, t) = c(k, t) / X^fcLi c (^> 



2. Traffic Models 



Our second example is in the context of traffic models. In this context, we discuss the so-called 
Bus Route Model (BRM) [20]. Here, one is interested in the initial conditions or parameters which 
result in a clustering of buses or a traffic jam. The model is defined on a 1 -dimensional lattice of 
size L. Each site i has two associated variables 7* and fa: (i) If a site i is occupied by a bus, = 1 ; 
otherwise n = 0. (ii) If a site i has passengers, then fa = 1; otherwise <^ = 0. A site cannot have 
both Ti = fa = 1, i.e., Tj + 0j < 1. If there are M buses, the bus density p = M/L is a conserved 
quantity. However, the total number of sites with passengers is not conserved. 

The update rules are as follows: (i) Pick a site % at random, (ii) If Tj = fa = 0, then set fa = 1 
with rate A, i.e., a passenger arrives at an empty site with rate A. (iii) If 7$ = 1 and r i+1 = 0, a bus 
hops onto a site with no passengers (fa+i = 0) with rate a, and to a site with passengers = 
1) with rate (3. Thus, the variables Tj — > and r i+1 — >■ 1 and — >■ with rate a or /3, as the case 
may be. Usually, < a as the buses slow down when passengers are being picked up. A jam in 
the system is a gap between buses of size x ~ O(L), which is stable in the thermodynamic limit. 

The MF approximation of this model considers the distribution of gaps P(x,t), ignoring the 
time-correlations in the hopping of buses. It should be noted here that, unlike the mass-transport 
models described in Sec. Ill AL the BRM is asymmetric. Thus, the movement of the buses is 
unidirectional although the hop rate is proportional to the size of the gap. These features put 
the BRM in a class of models which are referred to as zero-range processes (ZRP) 112 ill . The 
important property of a ZRP is that it yields a steady- state as a product of marginals calculated 
using well-defined procedures [22J]. Further, MF calculations are exact for this class of models. 

From simulations of the discrete model, heuristic arguments and MF theory, O'Loan et al. 
obtain evidence of a jamming transition as a function of the density of buses p. In terms of buses 
and passengers, the jam may be interpreted as follows. An ideal situation is one where the buses 
are evenly distributed along the route so that each bus picks up approximately the same number of 
passengers. Jamming or clustering of buses may occur if one of the buses gets delayed due to some 
fluctuation at a pick-up point. Subsequently, the buses which follow catch up with the delayed bus 
resulting in a jam! An important observation here is that the jamming is a consequence of a local, 
stochastic dynamics which couples the conserved variable (buses) and the nonconserved variable 
(passengers). The transition is reminiscent of the condensation transition described earlier ill ill , 
and has also been useful in describing clogging in the transport of sticky particles down a pipe 
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L20]. 



3. Granular Packing 

As a final example, we consider the packing of granular materials, which is important in many 
technological processes. The crucial issue in these problems is understanding the complex network 
of forces which is responsible for the static structure and properties of granular materials. One such 
system which has been subjected to experiments, simulations and analysis is a pack of spherical 
beads in a compression cell Q6Q . 

The bead pack is modeled as a regular lattice of sites, each having a particle of unit mass. The 
mechanisms which lead to the formation of force chains in the system are summarized in the rules 
defined below: (i) Each site i in layer D is connected to iV sites j in layer D + 1. (ii) Only vertical 
forces are considered explicitly. A fraction of the total weight supported by particle i in layer 
D is transmitted to particle j in layer D + 1. Thus, the weight W(D, i) supported by the particle 
at the i th site in layer D satisfies the stochastic equation 

W(D + = 1 + qa{D)W(D, i). (11) 

i 

The qij(D) are independently-distributed random variables which satisfy the constraint ^ ■ = 
1, required for enforcing the force-balance condition on each particle. 

In general, the values of W at neighboring sites in layer D are not independent. The MF 
approximation of this model ignores these correlations. Defining a normalized weight variable 
v = W/D, we want to obtain the force distribution Pd(v), i.e., the probability that a site at 
depth D is subject to a vertical force v. Within the MF approximation, it is possible to obtain a 
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found that, for almost all distributions of q, 



recursive equation for Pd{v). Coppersmith et al. 
the distribution of forces decays exponentially. However, a power-law decay was also observed in 
some cases. 

IE. FRAGMENTATION AND AGGREGATION OF fc-CHIPS 
A. 1-chip model 



Let us first consider the 1-chip model. The chipping kernel has the simple form 

g m (n) = w5 nA . 
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(12) 



With the above kernel, Eqs. CQ)-© become 

j m oo m2 

—P(m,t) = -P(m,t) J~] w8 mi>1 - P(m, t) ^ P(m 2 , t) ^ w5 mit 

7711=1 1712 = 1 7711 = 1 



+ P(m + mi, t)wS mi>1 + P(m -m 1 ,t) P(m 2 ,t)w5 mijl , m > 1, 

7711=1 mi=l 7712=77X1 

(13) 

, OO 7772 OO 

-P(0,t) = -P(0,t) ^ P(m 2 ,t) ^ ^,1+ P(m ls t)w(J roiil . (14) 

7772 = 1 7771 = 1 77ll = l 

Absorbing if into the definition of time, these equations simplify to the following form: 

^-P(m,t) = -(l + s 1 )P(m,t)+P(m + l,t) + siP(m-l,t), m > 1, (15) 
at 

^P(0,t) = -s x P(0,t) + P(l,t). (16) 

Here, we have defined si(t) = J2m=i P{ m i t) as tne probability of occupancy of a site with mass 
m > 1. Consequently, the probability of a site being empty is P(0, t) = 1 — si(t). 

The above rate equations were obtained earlier in Refs. [flll.[l2ll and were solved exactly. We 
recall this calculation to illustrate the generating-function approach for obtaining steady-state so- 
lutions of such rate equations. Defining the generating function Q(z, t) = J2m=i zm P(m, t), an 
equation for dQ/dt can be obtained from Eq. ([TBI by multiplying both sides by z m and summing 



over m: 



Pi P> 00 

-Q(z,t) = -$> m P(m,t) 



dt^ y ' ' dt 

777=1 



1 + si) ]T z m P{m, t) + J2 z m P{m + l,t) + Sl J2 z m P[m - 1, t) 

771=1 771=1 771=1 

1 OO OO 

1 + 8l )Q + z m P{m, t) + s x z z m P(m } t) 



z 

771=2 771=0 

= -(1 + Si)Q + -[Q- zP(l, t)) + Sl z [Q + P(0, t)) . (17) 

z 

Setting dQ/dt = 0, and substituting P(l) = si(l — si) from the steady-state version of Eq. (fT6l) . 
we obtain 

= £jilzLill£. (18) 

(1 - S\Z) 

The value of s% is fixed by mass conservation, which requires that J2m=i m P( m ) = P> where p is 
the mass density. Putting dQ/dz\ =1 = p, we obtain 



The steady-state distribution P{m) is the coefficient of z m in Q(z), and can be obtained by 
Taylor-expanding Q(z) about z — 0. This yields 



P(m) = (1 - s{)sf, m > 1. 



(20) 



For a more complicated function Q(z), we can obtain P{m) by inverting It is useful to 

illustrate this for the simple form of Q(z) in Eq. (IT8T ). Thus, 



P(m) 



dz 



Q{z) 



2ni J c " z m+1 ' 



m > 1. 



(21) 



Here, the closed contour C encircles the origin in the complex plane counter-clockwise and lies 
inside the circle | z \— l/s±. The integral is calculated using the residue theorem. Only those 
singular points which lie within C (viz., z = 0, which is a pole of order m) contribute to this 
evaluation. The associated residue is 



Res h{z = 0) 



1 



d 



m—l 



Q(z) 

(m - 1)! dz" 1 - 1 \ z 



2=0 



Thus, the steady-state mass distribution is 

p( m ) = -L.27riRes/i(0) = (l-si)s^, m >1. 

ZITl 

Notice that P(0) = 1 - s lf so Eq. dJU) is also valid for P(0). 
Using Eq. (fT9l ), the above mass distribution can be rewritten as 

1 



P(m) 



P 



1 + p \l + p 

ae- bm , 



(22) 
(23) 



(24) 



(25) 



where 



1 



In 



1 + P 



(26) 



1 + P V P 

In the case of simple chipping kernels, as in Eq. (fT2l) . the above solution can also be obtained 
directly from the difference equations (fT5l)-(fT6l) by setting the left-hand- side (LHS) to zero. We 
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can then write down expressions for the first few terms of P{m): 

P(l) = s x P(0) = si - si, 
P(2) = (l + Sl )P(l)- Sl P(0) 

= slP(0) = sl- si 
P(3) = (l + si)P(2)- Sl P(l) 

= slP(0) = st-s$, 

p{rn) = (1 + si)P(m - 1) - siP(m - 2) 

= s f P(0) = s™ - (27) 

which is identical to Eq. (|24l) . Again, the mass conservation condition Xlm=i m P{ m ) = P results 
in Eq. (fT9l) . as expected. 

The 1-chip solution in Eq. (1251) is important because of its universal nature. As a matter of fact, 
it is a steady-state solution for all MF models where the chipping kernel g m (n) is independent of 
the mass of the departure site m, g m (n) = g(n). To confirm this, we consider Eqs. CQ)-© with 
g m (n) replaced by g(n). The corresponding rate equations are 

, m oo ni2 

-P(m,t) = -P(m,t) g(mi) - P(m,t) ^ P(m 2 ,t) £ </(mi) 



mi=l m,2=l m\=l 

oo 



+ 2j + 



mi 



mi =1 



+ ^P(m-mi,t) P(m 2 ,t)g(mx), m > 1, (28) 



mi=l m2=m\ 



-P(0,t) = -P(0,i) ^ P(m 2 ,f) ^ffW + ]T PK,t) 9 K). (29) 

m,2=l mi=l mi = l 



In the steady state, the above equations may be combined to obtain 

m j x 

-P(m) ^ ^( mi ) _p( m )__ P(mi)g{mi) 

m\=l ^ ' mi=l 

00 00 00 

+ P(m + mi)g(mi) + P(m — mi ) ^ P(m 2 )s(mi) = 0. (30) 

mi=l mi=l m2=?7il 
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Substituting P{m) = aexp(-bm) on the RHS of Eq. (l30l) . we obtain 

m oo oo 

RHS = ]T - - J2 ae ~ bmi 9(™i) + £ ^^i) 

mi=l mi=l mi=l 

oo oo 

+ fle " t(m2 ^ l) 9( m i)' (31) 

mi=l m2=mi 

The first and fourth terms cancel, and the second and third terms cancel, so RHS = 0. This 
confirms that P{m) = ae~ bm is a solution of Eqs. (|28l) -(|29l). The constants a and b are fixed 
by the requirements of probability normalization Em=o P( m ) — 1] an d mass conservation 
Cm=i m -^'( m ) = Al- We leave it as an exercise to the reader to verify that these conditions 
lead to the same values of a and b as in Eq. (|26| ). 

Next, let us generalize the 1-chip model to a &;-chip model, where k > 1. These /c-chip models 
are interesting in physical situations where the deposited material consists of polymers or aggre- 
gates. We will see that the steady-state solutions for /c-chip models exhibit a k -branch structure. 



B. 2-chip model 

The steady-state distributions for the 2-chip model can be obtained using a procedure similar 
to the 1-chip model. The corresponding form of g m (n) is 

g m (n) = w5 ni2 - (32) 

The rate equations in this case are (absorbing w into time t) 

—P(m,t) = -(l + S2)P(m,t)+P(m + 2,t) + s 2 P(m-2,t), m > 2, (33) 
at 

^-P(m,t) = -s 2 P(m,t) + P(m + 2,t), m < 2. (34) 
at 

Here, s 2 (t) = Ylm=2^'( m ^) ^ s me probability of sites having mass 2 or more. Notice that the 
kernel in Eq. (l32l) is independent of the mass of the departure site. Thus, the 1-chip solution is 
a steady-state solution of Eqs. (I33l-d34l). as can be verified by direct substitution. However, an 
arbitrary initial condition P(m, 0) will not relax to this solution due to the presence of conserved 
quantities, as we shall see shortly. 

The steady- state generating function Q(z), obtained in analogy with the 1-chip case, is as 
follows: 

= sfr-^+^a-.,) (35) 

(1 - s 2 z 2 ) 
12 



where s 2 = J2m=2 P( m )- It is straightforward to Taylor-expand Q(z) in Eq. (1351) and identify 
P(m). Alternatively, we can obtain P{m) using Eq. (|2D) . Thus 

F(m) = / az ; rr , m > 1. (36) 

In this case, the contour C encircles the origin counterclockwise, and lies inside the circle | z \ = 
1 / y/s2. The singularities of the integrand f2(z) in the above equation are z = (pole of order m), 
z = l/y/~S2 (simple pole) and z = —1/y/s^ (simple pole). The second and third poles lie outside 
C, making Res ^(z = 0) the only contributing residue. This evaluation yields 

Res f 2 (z = 0) = 2-1^ [1 + (-l) m ] sf 2 + ^l^L [1 _ (_!)-] s { 2 m - l)/2 . (37) 

Thus, the steady-state probability distribution for the 2-chip model is given by 

P{m) = (1 - si)s^ /2 5 mod(m 2) + (si - s 2 )4 m ~ 1)/2( Wl( m ,2),i 

= P e (m) + P°(m). (38) 

Here the function mod(m, n) is defined as the remainder on division of m by n. The first and 
second terms on the RHS of Eq. (1381 ) are the steady-state distributions for even values of m [P e (m)} 
and odd values of m [P°(m)}, respectively. Thus, the 2-chip model has a steady- state solution 
comprising of two branches, both of which have the same exponential decay. Notice that the 
occupation probabilities for sites with even or odd units of mass are 

S e = Yl p6 ( m ) = T~~~^> ( 39 ) 

z — ' 1 — So 

m=0,2,4,.. z 

oo 

S a = pt » = T— ^- (40) 

m=l,i,b,.. 

These quantities remain conserved during the evolution because of the nature of the 2-chip move. 
The two branches appear as a consequence of these two conserved quantities. 

The probabilities of occupancy si and s 2 are related to the mass density p as follows: 



= dQ 

^ dz 



31 + 52 (41) 



*=1 1 - S 2 

As before, p may be calculated from either Q(z) or P(m). The quantities s\ and s 2 can be 
determined in terms of p and S s (or S ): 

P-S e + 1 p + S e -l 
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It should be noted that p + S e is always greater than I, ensuring s 2 > 0. If we choose the initial 
conditions P(m, 0) such that s 2 = sf in Eq. (l42l) . the branched solution of Eq. (1381) reduces to the 
1-chip exponential solution. Alternatively, if we substitute s 2 = si in Eq. (1351 ), we recover the 
generating function of the 1-chip model in Eq. (fT8l) . 



C. £>chip model 

In general, consider the case of k units of mass chipping from a site and then aggregating with 
the mass of a randomly-chosen nearest neighbor: 

g m (n) = w5 n)k . (43) 

The corresponding rate equations for P(m,t), obtained by substituting Eq. (|43l in Eqs. 
are as follows (absorbing w into t): 

^-P(m,t) = -(l + s k )P(m,t)+P(m + k,t) + s k P(m-k,t), m > k, (44) 
at 

^-P(m,t) = -s k P(m,t) + P(m + k,t), m<k. (45) 
at 

Here, s k (t) = Ylm=k P( m ,t) is the probability of sites having mass k or more. As the kernel 
in Eq. (|43l) is independent of the mass of the departure site, P{m) = aexp(—bm) is a steady- 
state solution of Eqs. (|44|)-(|451). However, as in the 2-chip case, an arbitrary initial condition 
P(m, 0) will not relax to this exponential solution due to the presence of k conserved quantities: 
^^ =0 P(n + mk, t) with n — 0, 1, k — 1. Rather, the steady-state solution will consist of k 
branches. 

The steady-state generating function for the /c-chip model is 

n(v \ g(fi ~ g 2) + z 2 (s 2 - s 3 ) + ■ ■ ■ ■ +z k s k (l - gi) 

Q{z) = ' m 

where s k = J2m=k P{ m )- The corresponding probability distribution is as follows: 

P(m) = (1 - si)s™ /fe <5 mod(mifc)i0 + (si - s 2 )s k m - 1)/k 5 mod{mm + .... + 

( Si - s l+ i)s k m - t)/k 5 mod{mtk)ji + .... + (s fc _i - Sfc)4 ro ~* +1)/ *^nod( ro ,*),fc-r ^ 

For completeness, we present the derivation of Eq. (|47l) for the 3-chip case in Appendix B. The 
Sj's are related to the mass density via the relation 

Si + So + + Sk 

= IT 2T (4g) 

1 - s k 
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Thus, the steady-state solution for the fc-chip model consists of k branches. All the branches in 
the probability distribution of Eq. (PTTT) decay exponentially with a slope \n(sk)/k. The quanti- 
ties si, , s k are determined from Eq. (1481 plus the (k — 1) conserved probability sums in the 

branches: 



m=i,i+k,... 
s i ~ s i+l 

1 - s fc 



0,1...., fc - 1. (49) 



Notice that one of the Si's (say, Sk-i) is not independent because Yli=o &i — 1 Further, appro- 
priately chosen initial conditions resulting in steady-state values s 2 = sf, s 3 = sf,....,Sfe = 
collapse the k branches in Eq. (|47|) to the 1-chip solution. 

Before concluding this subsection, let us make an observation about the A;-chip lattice model, 
i.e., the original model rather than its "rate equation" counterpart. This model has an exponen- 
tially large set of disjoint sectors, and configurations in different sectors are not connected by the 
dynamics. To see this, we denote the number of particles on a site j as m;. As k particles arrive 
at or leave this site at a given time, the quantity Mj = mod(mj, k) is conserved. Therefore, the 
set {Mi} is conserved by the dynamics, and labels a particular sector. The number of sectors is 
k N , where N is the number of lattice sites. Such systems have been referred to as many-sector 
decomposable systems by Menon et al. 12311 . 



IV. MONTE CARLO SIMULATIONS 



In this section, we present Monte Carlo (MC) results for some of the models discussed earlier. 
All simulations were performed on l-d and 2-d lattices with periodic boundary conditions. The 
lattice sizes were L = 1024 (in d = 1) and L 2 = 128 x 128 (in d = 2). The data presented here was 
obtained as an average over 500 independent runs. The details of the MC procedure are provided 
in Appendix C, so that the reader can implement these models numerically. 

First, we present results for chipping kernels which satisfy g m (n) = g(n), discussed at the end 
of Sec. lin Al In Fig. 2, we plot P(m) vs. m obtained from l-d MC simulations with three different 
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functional forms of g m (n): 



g m (n) = 1 



9m[n) = — 

-O.ln 



9 ' 

n z 



9m{n) = e- u - m (50) 

The MC data sets are numerically coincident with each other as well as the 1-chip solution in 
Eq. (|25l ) (denoted as a solid line), which was obtained from the corresponding MF equations. The 
mass density in each case was p = 5. There is excellent agreement between the different data 
sets, showing that the MC data is described very well by the solutions of the corresponding MF 
equations, even for d = 1. This feature is also observed in our subsequent results, suggesting that 
the MF equations are exact in the present context I 22L 

Next, we present results for the 2-chip model discussed in Sec. IIIIBI Figure 3(a) shows the 
steady-state distribution obtained from \-d and 2-d MC simulations for initial conditions with 
p = 10, and S & — 1, S — 0. The solid line denotes the result in Eq. (1381 ) with values of s 1 and 
s 2 evaluated from Eq. (l42l) . As our initial condition only had sites with even m populated, the 
steady-state solution is the even-m branch of Eq. (1381) . Figure 3(b) is similar, but for a mixed 
initial condition with p = 9.5 and S c = S Q = 1/2. The branched nature of the solution, resulting 
in a staircase-type probability distribution, is highlighted in the inset. 

In Fig. 4, we show the steady-state distributions obtained from l-d and 2-d MC simulations of 
the 3-chip model. In Fig. 4(a), the initial condition had P(9,0) = P(10,0) = P(11,0) = 1/3. 
This corresponds to p = 10, and all three branches are equally populated. The solid line in Fig. 4(a) 
denotes the result in Eq. (|47T) with s\, s 2 and s 3 calculated from Eq. (IB6I) . Figure 4(b) is analogous 
to Fig. 4(a), but the initial condition now has P(9, 0) = 1/2, P(10,0) = 1/3, and P(ll, 0) = 1/6. 
The corresponding value of the average density is p ~ 9.67. Again, in both sets of figures, the MC 
simulations agree very well with the corresponding MF result. 



V. SUMMARY AND DISCUSSION 



Let us conclude this paper with a summary and discussion. There has been much research in- 
terest in mass-transport models, which arise in many physical contexts. Therefore, we believe it is 
appropriate to make this subject accessible to a wider audience. This is the underlying motivation 
for this paper. For the purposes of this exposition, we focus on fragmentation-aggregation models 
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with conserved mass, i.e., there is no ongoing adsorption or desorption. We consider models with 
the chipping rate g m (n), where n is the chipped mass and m is the mass of the departure site. We 
use the corresponding mean-field (MF) equations to obtain the steady-state probability distribu- 
tions [P(m) vs. m] for different functional forms of g m (n). We show that a large class of chipping 
kernels, where g m (ri) is independent of m, give rise to an exponentially-decaying distribution: 
P{m) = aexp(-bm). This is also the MF solution for the 1-chip model [cf. Eq. (|25l) 1. where one 
unit of mass fragments from a site and aggregates with a randomly-chosen nearest-neighbor. 

We have also discussed A;-chip models for fragmentation and aggregation. The resulting steady- 
state distribution has k branches, each of which decays exponentially with the same slope. This 
slope is determined by the average density p, and the population of the branches in the initial 
condition P(m, 0) for the rate equations (Q])-©- The initial population in each of the branches is 
conserved during the evolution, and is also reflected in the steady-state distribution. 

Finally, we compared the MF analytical results with those from Monte Carlo (MC) simulations 
in d = 1, 2. In all cases, we found that the MF results for P{m) vs. m were in excellent agreement 
with the MC results. This demonstrates that the MF results are exact in the present context. 

There are many open research problems in the area of mass transport, aggregation and growth. 
These models can be tackled with a wide range of analytical and numerical techniques, which are 
both simple and elegant. We hope that this review will motivate further studies of this fascinating 
area. 
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Appendix A: Verification of the sum rule on P(m, t) 



We want to show that d Em=o P ( m ' *)] l dt = °- Usin § Ec l s - ©"©> 



d 
di 



£>(m,f) 



,m=0 



d 
It 



^P(m,t) + P(0,t) 



m=l 

oo 



11(2 



m oo 

= - P(m, t) Jm(mi) - *) 5^ P ( m 2, *) 9m 2 {^l) 

m=l mi=l m=l m2=l m\=l 

oo oo 

(mi) 

m=l mi=l 

oo m oo 

+ ^^P(m-m!,t) ^ P(m 2 ,t)g m2 (m 1 ) 

m=l mi=l rri2=mi 

oo rri2 oo 

-P(0,t) ^ P(m 2 ,t) 9m 2 (mi) + P(mx,t)g mi (mi)- ( A1 ) 

m,2=l m\=l mi=l 

We regroup terms and write 

oo m 

LHS = -^P(m,t) 9m{m x ) 

1=1 ?Tll = l 

OO OO 7T12 

- ^P(m,t) P(m 2 ,t) Y 9m 2 ( m i) 

m=l m2=l m\—\ 

oo r«2 

-P(0,0 ^ P(m 2 ,t) ^ 9mMx) 

m2=l mi=l 

oo oo 

22 + m i)%m+ mi (mi) + ^ P(m 1 ,t)g mi (m 1 ) 

.m=l mi=l mi=l 

oo m oo 

+ J2 J2 P(™~™i,t) Yl P(m2,t)g m2 (m 1 ) (A2) 

ra=l m\=l m,2=m\ 
oo m oo oo m,2 

- Y P(m, t) ^2 ft»W - P ( m ' p ( m2 ' *) $Z 9m 2 (mi) 

m=l mi=l m=0 rri2=l mi=l 

oo oo 

+ ^ 2J + m l^)^"i+mi(^l) 

m=0 mi=l 
00 m 00 

+ ^^P(m-mi,t) ^ P(m 2 ,t)^ m2 (mi). (A3) 



+ 



+ 



m=l mi=l 



rri2=rn\ 



Our initial condition for Eqs. ©-© satisfies Ylm=o P{ m i 0) = 1- Therefore, we set 
J2m=o P( m i i) = 1 on the RHS of Eq. (|A3b . This is justified subsequently as it results in 



18 



d\E™ =Q P(m,t)]/dt = 0. Thus 

oo m oo m,2 

g m2 (mi) 

m=l mi=l 7712=1 mi=l 
oo oo 

m2\ m l) 

mi = l 7712 = 7711 

OO OO OO 

+ 5^ ^ P(m-m 1; t) P{m 2 ,t)g m2 {m x ) (A4) 

mi = l 771=7771 7772=7771 

OO 777 OO OO 

= -2j2 p (m,t) 9m{m x )+2 J2 P{™2,t)g m2 (m x ). (A5) 

777=1 7711 = 1 7771=1 7712=7711 

In the second term on the RHS of Eq. (IA5I) . we interchange the order of the summations over m x 
and m 2 . This leads to a cancellation of the two terms, proving that d Em=o P( m i t)) /dt = 0. 
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Appendix B: Steady-State Distribution for the 3-chip Model 



The generating function for the 3-chip model is [cf. Eq. ([461) 1 

= fr-sjz+isi-s^ + ssil-sjz* (B1) 
(1 - S 3 Z d ) 

where s 3 = J2m=3 P( m )- From Eq. (1271) . we have 

vt \ 1 f a ( Sl ~ + ( g2 ~ s ^ z<2 ± ^ ~ Sl )^ 3 >> 1 

_r m = / dz — -. , m > 1. (B2) 

V 7 2vri J c z m+1 {l-s 3 z 3 ) ' V 7 

The closed contour C encircles the origin and lies inside the region defined by | z |< I/S3 1 / 3 . As 
usual, only Res f 3 (z = 0) contributes to the integral, and is evaluated as 



Res f 3 {z = 0) 



d 



m—l 



Q(z) 

(m - 1)! dz" 1 - 1 \ z 



z=Q 



1 d 



m—l 



(m - 1)! dz" 1 - 1 



(si - s 2 ) + (s 2 - s 3 )2: + s 3 (l - si)z 2 



-s 3 ) (z* - l/s 3 ) 



z=0 



.3 

(fl - gg) 

3 

(«2 - S3) 



1 + 2 cos 



2?r(m - 3) 



^m/3 



1 + 2 cos 



1 + 2 cos 



2vr(m - 1) 



2?r(m - 2) 



,(m-l)/3 



Jm-2)/3 



(B3) 



Then, the steady-state mass distribution is 



P(m) 



(l-*i) 
3 



1 + 2 cos 



27r(m - 3) 



^m/3 



(>2 - S3) 



1 + 2 cos 



1 + 2 cos 



27r(m - 1) 



27r(m - 2) 



(m-l)/3 



(m-2)/3 



(1 - Si)s 3 m/3 <J m od(m,3) 1 + («1 - S 2 )s 3 (m 1)/3 5mod(m,3), 
+ (S 2 - S3)s3 {m " 2)/3 5 mo d( m ,3),2- 



(B4) 



The relation between the mass density p and the Sj's can be obtained directly by using Eq. (|B4I) . 
Thus 



Si + s 2 + g 3 



(B5) 
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As for the 2-chip model, the s/s can be obtained as a function of p and the probability sums on 
each of the branches in Eq. (|49| ). These result in the following expressions: 

2S\ + S0S2 — S±S 2 



si = . n ' „ ~ + 

S2 



p + 2 + S -S 2 (p + 2 + S - S 2 )(l - S 2 ) 

p S 2 - Si 

p + 2 + S -S 2 p + 2 + S , -S 2 ' 



p Si + S 2 + SqS 2 — S 2 . , 

H ~ p + 2 + S -S 2 ~ (p + 2 + S -S 2 )(l-S 2 )' { ' 
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Appendix C: Monte Carlo Simulations of Mass-Transport Models 
1. 1-Chip Model 

Consider a 1 -dimensional lattice with periodic boundary conditions and label the sites with 
integers i = 1,2, 3, .., L. Fix the mass density p. 

1 . Initializing the lattice: Integer masses are placed on the lattice sites in accordance with 
the chosen p such that J2i=i m i = PL- A simple procedure for achieving this is as follows: 

(a) Choose an integer random number NRAN(i) from the range [0, L] for each site i. 

(b) Assign m, = Int [{p + 1)L * NRAN(?)/ £\ NRAN(i)], where lnt(x) refers to the in- 
teger part of x. 

2. Chipping and aggregation: A site % is chosen at random. If rrii is non-zero, a unit mass chips 
and aggregates with a randomly-chosen neighbor. Thus mj — >■ to* — 1, and mj_i — > + 1 
or m i+ i — > m i+ i + 1 with a probability 1/2. 

3. Repeat step 2 for L times, which corresponds to one MC step (MCS). 

4. Compute the mass distribution of lattice sites. 

5. Repeat steps 2-4 for several MCS, storing the mass distribution at intermediate MCS. 

6. Repeat steps 2-5 for a large number of independent lattice configurations, generated via step 
1. 

7. Compute the configuration-averaged steady-state mass distribution, P(m) vs. m. 

The steady state is achieved when the mass distribution does not change (apart from numerical 
fluctuations) in subsequent MCS. 

2. Case with g m (n) = 1/n 2 

In general, the fragmentation rates g m (n) can be computed prior to the simulation and stored 
in a matrix for easy look-up using the following tips: 
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1. Define a matrix G of dimension N x N, chosen specific to g m (n). For example, if g m (n) = 
1/n 2 , the rate of chipping a mass of size 1000 is 1CT 6 . We can assume that chipping of 
masses greater than 1000 occurs with a very small rate, and hence these events may be 
ignored. Thus we may set N = 1000. 

2. Initialize G mn to zero. The row index m corresponds to the mass of the lattice site chosen to 
fragment. Masses greater than N may be treated as iV for the computation of fragmentation 
rates for reasons discussed in 1 above. The column index refers to the chipped mass n < m. 
Thus the matrix G has a triangular form, with the non-zero entries calculated using g m (n). 
For example, if g m {n) = 1/n 2 , the rates for the 10 th row are [1, 0.25, 0.1111, 0.0625, 0.04, 
0.02777, 0.0204, 0.0156, 0.0123, 0.01, 0, 0, 0]. To connect these to probabilities we 
normalize these numbers by Yl^Li V™ 2 - 

Once the look-up table for fragmentation rates is computed, the MC procedure is as before with 
step 2 replaced by the following steps. 

A. A site % is chosen at random. If rrii is non-zero, draw a random number r in the interval 
(0,1). 

B. Go to the mf 1 row of the table and check the two consecutive entries which sandwich r. 
The column number of the larger entry is the number of mass particles n that chip from m*. 
Thus Tfi-i — y 777-j — n, an d mj-i — > m^x +nor m i+1 — > m i+1 + n with probability 1/2. 
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FIG. 1: Schematic representation of the conserved-mass model with fragmentation and aggregation. A 
mass n can chip from a site with mass m with a fragmentation rate g m {n), and aggregate with the right 
(left) neighbor with probability p (1 — p). 
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FIG. 2: Steady-state probability distributions [P(m) vs. m] from l-d MC simulations with three different 
functional forms for the chipping kernel g m (n). The data sets are plotted on a linear-logarithmic scale. The 
details of the MC simulations are provided in the text. The mass density for the initial conditions is p = 5. 
The solid line denotes the 1-chip solution in Eq. (|25T ). 
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FIG. 3: Steady-state distributions for the 2-chip model, obtained from MC simulations, (a) Plot of P(m) 
vs. m on a linear-log scale, from MC simulations in d = 1,2. The initial conditions were characterized 
by average density p = 10, and S e = 1, S Q = 0. The solid line denotes the solution in Eq. (I38T ) with 
si = S2 = p/(p + 2). (b) Analogous to (a), but the initial conditions for the MC simulations have a mixture 
of both even and odd masses. The corresponding parameter values are p = 9.5, S c = S = 1/2. The solid 
line denotes the solution in Eq. OH) with s x = (2p + l)/(2p + 3), s 2 = (2p - l)/(2p + 3). The inset 
highlights the staircase structure of the probability distribution. 



28 




FIG. 4: Analogous to Fig. 3, but for the 3-chip model, (a) Plot of P(m) vs. m from MC simulations with 
P(9, 0) = P(10, 0) = Pill, 0) = 1/3, so that p = 10. The solid line denotes the result in Eq. §7) with 
s\, S2 and S3 calculated from Eq. (IB6b as si = 0.917, S2 = 0.833, S3 = 0.75. (b) Plot of P{m) vs. m for 
initial conditions with P(9,0) = 1/2,P(10,0) = 1/3,P(11,0) = 1/6, so that p ~ 9.67. The solid line 
denotes the result in Eq. (@7]> with si = 0.875, s 2 = 0.792, s 3 = 0.75. 
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